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Abstract 

Planetary topography can either be modeled as a load supported by the lithosphere, or as a 
dynamical effect due to lithospheric flexure caused by mantle convection. In both cases the response 
of the lithosphere to external forces can be calculated with the theory of thin elastic plates or shells. 
On one-plate planets the spherical geometry of the lithospheric shell plays an important role in the 
flexure mechanism. So far the equations governing the deformations and stresses of a spherical shell 
have only been derived under the assumption of a shell of constant thickness. However local studies 
of gravity and topography data suggest large variations in the thickness of the lithosphere. In this 
article we obtain the scalar flexure equations governing the deformations of a thin spherical shell with 
variable thickness or variable Young's modulus. The resulting equations can be solved in succession, 
except for a system of two simultaneous equations, the solutions of which are the transverse deflection 
and an associated stress function. In order to include bottom loading generated by mantle convection, 
we extend the method of stress functions to include loads with a toroidal tangential component. We 
further show that toroidal tangential displacement always occurs if the shell thickness varies, even 
in the absence of toroidal loads. We finally prove that the degree-one harmonic components of 
the transverse deflection and of the toroidal tangential displacement are independent of the elastic 
properties of the shell and are associated with translational and rotational freedom. While being 
constrained by the static assumption, degree-one loads can deform the shell and generate stresses. 
The flexure equations for a shell of variable thickness are useful not only for the prediction of the 
gravity signal in local admittance studies, but also for the construction of stress maps in tectonic 
analysis. 

1 Introduction 

Terrestrial planets are flattened spheres only at first sight: close-up views reveal a rich topography with 
unique characteristics for each planet. Unity in diversity is found by studying the support mechanism 
for topographic deviations from the hydrostatic planetary shape. A simple mechanism, called isostasy, 
postulates mountains floating with iceberg-like roots in the higher density mantle. Another simple mech- 
anism assumes that mountains stand on a rigid membrane called the mechanical lithosphere. The two 
simple models predict very different gravity signals, from very weak in the first to very strong in the 
second. The truth lies in-between: an encompassing model views topographic structures as loads on an 
elastic shell with finite rigidity, called the elastic lithosphere (the elastic lithosphere is a subset of the 
mechanical lithosphere). The rigidity depends on the elastic properties of the rocks and on the apparent 
elastic thickness of the lithosphere. The latter parameter is the main objective of many studies, since its 
value can be related to the lithospheric composition and temperature. Another important application of 
the model of lithospheric flexure is the determination of stress maps, which can then be compared with 
the observed distribution of tectonic features. 
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Besides the important assumption of elasticity, the model of lithospheric flexure is often simplified 
by two approximations. The first one states that the area to be analyzed is sufficiently small so that the 
curved lithospheric shell can be modeled as a flat plate. The second approximation states that the shell 
(or plate) is thin, which means that deformations are small and that the elastic thickness is small with 
respect to the wavelength of the load. 

On Earth, the model of lithospheric flexure has been very successful for the understanding of the 
topography of the oceanic plates, whereas the analysis of continental plates is fraught with difficulties 
due to the very old and complex structure of the continents. As far as we know, plate tecton ics do not 
occur at the present time on other terrestrial planets, which are deemed one-plate planets 



Solomoi 



1978J , although the term single-shell would fit better because of the curvature. A single shell can support 
loads of much larger extent and greater weight, the best example of which is the huge Tharsis volcanic 
formation covering a large portion of Mars. Indeed such a load cannot be supported by the bending 
moments present in a thin flat plate, whereas it can be supported by stresses tangent to the shell: the 
shell acts as a membrane. 

1969| for 



The first application of thin shell theory to a pl anet was done by 



Brotchie and Silvester 



the lithosphere of the Earth (see also 



Brotchie 



197l|). However the approximation of fiat plate theory 
was se en to be sufficien t when it was understood that the lithosphere of the Earth is broken into several 
plates 



Tanimoto 



Brotchie and Silvester 



19691 ] do not consider tangential loads and their flexure 
eq uation only in clud es dominant te rms in derivatives; their equation is thus a special case of the equations 
of 

theory of Kraus or Vlasov unless ment i oned otherwise. 

Solomon and Head ation to study displacement and stress 



Kraus 



19671 and 



Vlasov 



1964j discussed below. The articles reviewed hereafter use the thin shell 



On the Moon, 



Turcotte et al. 



in mare basins, 
the type of stress suppor ting topographic loads. 



198l| estimated g ravity-topography rat ios for the mascons and discussed 



Arkani-Hamed 



1998| mode l ed th e support of mascons 



with Brotchie's equation. \Suaano and Hekn 2004] and I Crosby and McKenzieX 20051 ] estimated the elastic 
thickness of the lithosphere from Lunar Prospector data. 

On Mars, the dominance of the Tharsis rise in t he topography led to num e rous applicat i ons o f 
the theory of thi n ela stic shells to lithospheric fl exure. 



Thurber and Toksdz 



Hall et al 



19861 and 



Janle and Jannsen 



thickness under Martian volcanoes. 
memb rane regimes. 



1978 1 . 



Comer et al 



19861 us e d Bro tchie's equation to estimate the lithospheric 



Turcotte et al 



Willemann and Turcotte 



19811 ] studied the transition between ben ding and 



1982j analyzed the litho spheric support of Tharsis. \Sleev an d Phillips 



1985] ana lyzed the membrane stress distribution on the whole surface. 



Banerdt et al 



1992J, 



Banerdt and Golombek 



2000| | and lPhillivs et ali 2001 1 used a mo del of lithospheric flexure including membrane s upport, bending 



stress es and tangent loads 



Banerdt 



1986| in order to study the global stress distribution. 



Arkani-Hamed 



Johnson et al 



2000] determined the elastic thickness beneath large volcanoes with Brotchie's equation whereas . 
20001 estimated the ela stic thickness beneath th e North Polar Cap. Using local admittance analysis with 
spat iospectral methods. \McGovern et all 1200211 det ermined the elastic thickness at various locations [see 



McGovern et al 



2004 1 . 



McKenzie et al 



20021 ] made local admittance ana l yses o f line-of-sight grav 



also 

ity data both with flat plates and with spherical shell models. 

shell formula for a one-dimensional wavelet analysis of the adm itta nce in order to determine the average 



Turcotte et al 



200 



3] used the spherical 



elastic thickness of the lithosphere. \Zhona and Roberts 



200.^ 



and 



Lowru and Zhona 



2003] studied the 



support of the Tharsis rise with an hybrid model i ncluding the flexure of a thin elastic shell as well as the 
internal loading of a thermal plume in the mantle. Bellequic et al. 2005] determined the elastic thickness 
and the density beneath large volcanoes. 



Searls et al. 



20061 ] investigated the elastic thickness and the 



density beneath the Utopia and Hellas basins. 

Venus is considered as a one-plate planet but does not have giant volcanic or tectonic structures 
comparable to Tharsis. The spherical shell model has thus not been used as often for Venus as for Mars. 
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Banerdt 



1986] studied the global stress distribution. 



Janle and Jannsen 



1988j | and 



Johnson and Sandwell 



1994 | used Brotchie's equation to estimate the lithospheric thickness in various lo cations. 



Sandwell et al 



1997] computed global strain trajectories for comparison with observed tectonics. \Lawrence and Phillivs 



20031 ] inverted the admittance in order to esti mate the elastic thickness and mantle density anomalies 



Anderson and Smrekar 



200 



(J established a global map 



over two lowland regions and one volcanic rise. 

of the elastic thickness based on local admittance analysis. Mercury's topography and gravit y fields are 
not ye t known well enough to warrant the application of a thin elastic shell model. We refer to 



Wieczorek 



20071 ] for a review of recent results regarding the lithosphere of terrestrial planets. 



The elastic thickness o f the lithosphere is no t at all homogeneous over the surface of a planet. For 



McGovern et al 



2004] (for Mars) and \Anderson and Smrekan [20061 ] (for Venus) find litho 



example, 

spheric thickness variations of more than 100 km. The former study explains the variation in lithospheric 
thickness in terms of different epochs of loading, as the lithosphere is thickening with time. Other studies 
however in ferred that spatial variations in lithospheric thickne ss on Mars are as important as temporal 
variations 



Solomon and Head , 11982 , 



1990; 



Comer et al. 



1985]. Thin spherical shell models have always 



been applied with a constant elastic thickness for the whole lithosphere. Local studies ar e done by win 
dowing both the data (gravi t y and topography) and the model predictions for gravity 



1997; 



Wieczorek and Simons 



Simons et al 



20051 ]. The assumptions behind these methods are that the elastic thick- 



ness is constant within the window and that the area outside the window can be neglected. The first 
assumption is of course true for a small enough window, but the size of the window is limited from below 
by the resolution of the data 



Wieczorek and Simons 



2005]. Even if the first assumption were true, the 



second assumption is violated in two ways (unless the elastic thickness is spatially constant). First, the 
deformation of the shell within the window as well as the associated stress field are both modified if the 
elastic thickness is changed in the area outside the window. Second, the value of the predicted gravity 
field within the window depends on the shell deflection outside the window. 

These reasons make it interesting to develop a model of the lithospheric flexure for a spherical shell 
of variable thickness. Although a full inversion of the gravity and topography data is impractical with 
such a model because of the huge size of the parameter space, other applications are of high interest. For 
example a two-stage inversion can be considered: in the first stage a constant elastic thickness is assumed, 
and the resulting values are used in the second stage as a starting point for an inversion with variable 
elastic thickness (the parameter space can also be constrained with an a priori). Moreover this model 
can be used to produce synthetic data and thus allows us to check the validity of inversions assuming a 
constant elastic thickness. Finally, stress and strain fields can be computed for given variations of the 
elastic thickness, with the aim of comparing stress and strain maps with tectonic features. 

General equa t ions governing the deformations of a thin elastic shell have been given by various 
authors [e.g. 



Love 



1944; 



Vlasc 



1964; 



Kraus 



19671 ] . However the possibility of variable shell thickness 



is only considered at the early stage where the strain-displacement relationships, Hooke's law and the 
equilibrium equations are separately derived for the thin shell. The combination of these three sets of 
equations into a unique equation for the transverse deflection is made under the restrictio n of constant 
elastic thickness, 'owing to the analytical complications which would otherwise arise' [see I Krausv 11967 . 
p. 199]. 

This article is dedicated to the derivation of the minimum set of equations governing the deformations 
of a thin elastic spherical shell with variable thickness. Using Kraus' method of stress functions, we find 
that the transverse deflection is the solution of a simultaneous system of two differential equations of the 
fourth order. Contrary to the case of constant thickness, these equations cannot be combined due to 
the presence of products of derivatives of the thickness and derivatives of the deflection or of the stress 
function. In order to include bottom loading generated by mantle convection, we e xtend the me thod of 
stress functions to include toroidal tangential loads, which were not considered by I Krausl 19671 ]. Non- 



3 



toroidal loading is for example generated by tangential lithostatic forces whereas toroidal loading could 
be due to mantle flow generating drag at the base of the lithosphere (mantle flow also produces non- 
toroidal loading) . With applications to tectonics in mind, we derive the equations relating the tangential 
displacements and the stresses to the transverse deflection and the stress functions. We further show 
that toroidal tangential displacement occurs even if there is no toroidal loading (unless the shell thickness 
is constant). Finally we prove three properties specific to the degree-one harmonic components: (f) 
the degree-one transverse deflection and the degree-one toroidal tangential displacement drop from the 
elasticity equations because they represent rigid displacements of the whole shell, (2) the transverse and 
tangential components of degree-one loads are related so that the shell does not accelerate, (3) degree- 
one loads can deform the shell and generate stresses. Though our aim is to introduce a variable shell 
thickness, the final equations are also valid for a variable Young's modulus (Poisson's ratio must be kept 
constant). 

Another way to take into account variations of the lithos pheric thickness cons ists in treating the 
lithosphere as a th r ee-dimensional s pheri c al solid that is elastic [Metivier et al. 



Zhona et al. 



2003; 



Latvchev et al 



2005: 



20061 ] or viscoelastic [e.g. 



Wang and Ww , l2006j . The resulting equations are exact (no 



thin shell approximation) and can be solved with finite clement methods. The thin shell assumption 
is probably satisfied for known planetary lithospheres; in case of doubt, it is advisable to compare the 
results of thin shell t heory with three-dimensional models assuming constant elastic thickness: static 
thick shell models [e.g. 



Banerdt et al 



f982: 



Janes and Melosh 



f990: 



2000; 



Arkani-Hamed and Riendler , 



Zhona 



2002 



2002J. The advantage of thin 



or time-dependent viscoelastic models [e.g. \Zhona and Zuben 
shell equations is their two-dimensional character, making them much easier to program and quicker to 
solve on a computer. Solving faster either gives access to finer two-dimensional grids or allows to examine 
a larger parameter space. 

We choose to work within the formalism of differential calculus on curved surfaces, without which 
the final equations would be cumbersome. Actually the only tool used in this formalism is the covariant 
derivative, which can be seen by geophysicists as just a way of combining several terms into one 'deriva- 
tive'. All necessary formulas are given in the Appendix. The possibility of writing a differential equation 
in terms of covariant derivatives (in a tensorial form) also provides a consistency check. The presence of 
derivatives that cannot be included into covariant derivatives is simply forbidden. This is not without 
meaning for differential equations of the fourth order including products of derivatives. 

In section [2j we show how to obtain the strain-displacement relationships, Hooke's law and the 
equilibrium equations for a thi n spherical s hell. These equations are available in the literature for the 



general case of a thin shell [e.g. \Krausl Il967| . but we derive them for the spherical case in a simpler way, 
starting directly with the metric for the spherical shell. We examine in detail the various approximations 
made to obtain the thin shell theory of flexure, refraining until the end from taking the 'thin shell' 
quantitative limit in order to ascertain its influence on the final equations. In section [3l we use the 
method of stress functions to obtain the flexure equations governing the displacements and the stresses. 
In section^ we give the final form of the flexure equations in the thin shell approximation. We also study 
the covariance and the degree-one projection of the flexure equations. In section [5j we examine various 
limit cases in which the flexure equations take a simpler form: the membrane limit, the Euclidean limit 
and the limit of constant thickness. 
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2 Fundamental equations of elasticity 



2.1 Three-dimensional elasticity theory 

Linear elasticity theory is based on three sets of equations. We directly state them in ten sorial form for 



Ranalli, 



1987 



an isotropic mate r ial, si nce they are derived in Cartesian coordinates in many books [e.g 
Svnae and Schilm . Il978| . Recall that, in tensorial notation, there is an implicit summation on indices 
that are repeated on the same side of an equation. The first set of equations includes strain-displacement 
relationships: 



where is the infinitesimal strain tensor and Ui are the finite displacements. The 'comma' notation 
denotes the spatial derivative (see Appendix 17. 2j) . 

The second set includes the constitutive equations of elasticity or Hooke's law, relating the strain 
tensor and the stress tensor <7;,-: 



XeSa +2G 6A 



(2) 



where e = en + £22 + £33 and #y = 1 if i=j, otherwise it equals zero. The parameter A is known as the 
first Lame constant. The parameter G is known as the second Lame constant, or the shear modulus, or 
the modulus of rigidity. Boundary conditions are given by 



T, 



(3) 



with nj being the normal unit vector of the surface element and Ti being the surface force per unit area. 

The third set includes equations of motion which reduce to equilibrium equations for stresses if the 
problem is static: 

= • (4) 

Body forces, such as gravity, are assumed to be absent. Both strain and stress tensors are symmetric: 
€ij — tji and o~ij = o~ji. 

These three-dimensional equations do not yet have the right form for the description of the defor- 
mations of a two-dime nsio nal spherical shell. Various methods have be en used to generate appropriate 
equations. 



Love 



1944| and I Timoshenko and Woinowsky-Krieaen 1964j derive strain-displacements and 



equilibrium equations directly on th e surface of t he sp here (the latter only for the special cases of no 
bending or axisymmetrical loading). \Sokolnikoff\ 10561 ] derives strain-displacement equations in three- 
dimensional curvilinear coordinates using an arbitrary dia gonal metric but states without proof the equi- 
librium equations in curvilinear coordinates. \Kraus\ 19671 ] uses Sokolnikoff's form of strain-displacement 



equations and derives equilibrium equations for an arbitrary two-dimensional surface using Hamilton's 
principle (i.e. virtual displacements). 

Instead of directly deriving equations on the two-dimensional surface of the sphere, we will first 
obtain their form in three-dimensional curvilinear coordinates and then restrict them to the surface of 



the sphere. The first step can elegantly be done through the use of tensors \Synae and Schildx . I197S 
Equations ((T|)- (J4j) are tensorial with respect to orthogonal transformations, but not with respect to other 
coordinate transforms (one reason being the presence of usual derivatives). In other words, they are 
only valid in Cartesian coordinates. In a three-dimensional Euclidean space, tensorial equations have 
a simplified form in Cartesian coordinates because supplementary terms that make them tensorial with 
respect to arbitrary coordinate transformations are zero. The missing terms can be reconstructed by using 
a set of rules, such as the replacement of usual derivatives by covariant derivatives and the substitution 
of tensorial contraction to sum on components. Correspondence rules lead to the following three sets of 
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equations: 

eij = i (uiu + u^i) , (5) 

a tj = Xeg tj + 2Ge ij , (6) 

9 3k ^j\k = 0, (7) 

where e = g kl eki- The notation Uju denotes the covariant derivative of Ui (see Appendix 17. 2| . The 
metric and its inverse are noted gij and g lJ , respectively. Tensorial components cannot be expressed in 
a normalized basis (except for Cartesian coordinates) which is more common for physical interpretation 
(see Appendix 17. ip . Covariant components in equations ([5])- ([7]) are related to components defined in a 
normalized basis (written with a hat) by: 



ga Ui , 



\J 9ii9jj e ij 



&ij '\Z9ii9jj &ij i 

where there is no implicit summation on repeated indices. 

In the next section we will introduce additional assumptions in order to restrict the equations to the 
two-dimensional surface of a spherical shell. 

2.2 Spherical shell 

2.2.1 Assumptions of the thin shell theory 

Suppose that the two first coordinates are the colatitude 9 and longitude ip on the surface of the sphere, 
whereas the third coordinate £ is radial. R is the shell radius. 



Assumptions of the thin shell theory are [see I Kraud . 119671 chap. 2.2]: 

1. The shell is thin (say less than one tenth of the radius of the sphere). 

2. The deflections of the shell are small. 

3. The transverse normal stress is negligible: = 0. 

4. Normals to the reference surface of the shell remain normal to it and undergo no change of length 
during deformation: cqq = = = 0. 

The second assumption allows us to use linear equ ations to desc ribe the deflections. The third and fourth 



assumptions are not fully consistent: we refer to I Kraud 1967] for more details. We will relax them in 



the derivation of the equations for the deflection of a spherical shell. The crucial assumption is cr^ = 
which is essential for the restriction of Hooke's law to the two-dimensional shell. As we will see later, a^^ 
cannot be zero since it is related to the non-zero transverse load (besides the fact that it is incompatible 
with a vanishing transverse strain). What is absolutely necessary is that cr^ <C era for i = (0,<p). In 
section 15.3. 3[ we will show that this condition is satisfied if the wavelength of the load is much larger 
than the thickness of the shell. 

The reference surface is the middle surface of the shell. With the aim of integrating out the third 
coordinate, a coordinate system is chosen so that the radial coordinate £ is zero on the reference surface. 
The metric is given by 

ds 2 = (R + Cf (d9 2 + sin 2 dip 2 ) + dC, 2 . (8) 
Christoffel symbols necessary for the computation of the covariant derivatives are given in Appendix l7.3l 



(i 



2 R + 


( 


a C,C» 




1 1 




2 R + 


C 


1 1 




2 R + 


C 



2.2.2 Strain-displacement equations 

With the metric |8]), the strain-displacement equations ([5]) become 

1 

eee = p , ~ t (ug.g + u c ) > 
-ft + C, 

e w = 75 t (csc0Um„ + cot#u + fi^) 

H + (, 

ee v = - p | t (csc 8 ue,<p — cot8u v + u {p ,g) , (9) 

e^C = ((i? + C)^,C~^ + csc6l?i C^) 

where all quantities are given in a normalized basis. The fourth assumption of the thin shell theory 
implies that the displacements are linearly distributed across the thickness of the shell, with the transverse 
displacement being constant. Displacements can thus be expanded to first order in £: 

(u$, u Vl U() = {ve + C/3e,v v + C,(3 V , w) . (10) 

The coefficients (vg, v v , w, f3g,fi<p) are independent of C : (ve,v v ,w) represent the components of the 
displacement vector of a point on the reference surface, whereas (f3g, (3 V ) represent the rotations of tangents 
to the reference surface oriented along the tangent axes. We determine (3g and (3 V by applying the fourth 
assumption of the thin shell theory and the expansion (|10[) to the last two strain-displacement equations 
©: 

Pe = (ve - wfi) , 

P v = —(vtp-cscOw^). (11) 
R 

The substitution of the expansion (|10p into the fourth strain-displacement equation (e^ = uqx) l ea ds to 
the condition = 0, as postulated by the thin shell theory. 

After the insertion of the expansion formulas (fT0|) and (fTTj) , the first three strain-displacement equa- 
tions © become 

Zee = ^o + YTclR 4 ' 



C 



2e ev = l°e v + 1 + c/R 
The extensional strains ej), e° and jg are defined by 

e e = ft (v$,e + w) , 



e° = - 



(csc 8 v ViV + cot 9 Vg + w) , (13) 



lev = ~B ( v vfi - cot ^ v <p + csc 9 v e, v ) 



R 



They represent the normal and shearing strains of the reference surface. The flexural strains Kg, k£ and 
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t° are given by 



1 

Ii 2 
1 



-152 °i w > 



(14) 



They represe nt the changes in curvature and the torsion of the reference surface during deformation 



Kraus 



1967]. The differential operators ©1,2,3 are defined by 



dd 2 4 
csc 2 



O-x = esc ( 



1, 

d 2 
dip 2 
( d 2 



cote 



d_ 



1 



(15) 



— cot 1 



d_ 



The zero upper index in (ejj,e° 7$ , k 



K o _o 



\d6dip 

refers to the reference surface and appears in order to 



follow Kraus' notation. 

In Love's theory, one makes the approximation = in equations (|12p . However this does not 
simplify the calculations when the shell is a sphere because it has only one radius of curvature (for a 
shell with two radii of curvature, this approximation is a great simplification). Our choic e to keep the 
factor le ads to the same re sults as the theory of Fliigge-Lur'e-Byrne, explained in I Kraud [1967, 
chap. 3.3a] or \Novozhiloi\ 1964 , p. 53], in which this factor is expanded up to second order. In any 



case the choice between the approximation or the expansion of this factor does not affect the equations 
(derived in the next section) relating the stress and moment resultants to the strains: they are the same 



if 



-R+C 



is approximated to zeroth order, expanded to second order or fully kept. 



2.2.3 Hooke's law 

When the metric is diagonal and the basis is normalized, Hooke's law ([6]) becomes 



(Mi). 



fe=i 



There is no implicit summation on repeated indices. 

The third assumption of the thin shell theory, = 0, can be used to eliminate from Hooke's 

law: 

E 



l-I/ 2 

E 



{egg + ve w ) , 
2 (<W + v ^se) ■ 



1 — V 

Young's modulus E and Poisson's ratio v are related to Lame parameters by 

Ev 

A 
G 



(1 + i/)(1-2i/) ' 
E 



2(l + i/) 



In principle, the fourth assumption of the thin shell theory, eg^ = = 0, leads to ag^ — <t v q = but 
non-vanishing values must be retained for purposes of equilibrium. 
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M cpe 

Figure 1: Stress resultants and stress couples acting on a small element of the shell. The directions of the 
stress resultants (simple arrows) and the rotation sense of the stress couples (double arrows) correspond 
to positive components (tensile stress is positive). Loads (qe,qip,o) act on the reference surface. 



The substitution of the expansion (|12p into the thin shell approximation of Hooke's law gives 

(«8 + *"$)) , 



E ( o o C 

aee = v e + ve ^ + T+c/R 



E + + ^77* (< + »$)) > (16) 



i - 1/ 2 V v l + C/-R ' 

g ^0 | C _Q 



with eg, e°, 7$ , Kg, k° and r° defined by equations (JT3J) and (fl4|). The expressions for &g^ and a V Q will 
not be needed. 

We now integrate the stress distributions across the thickness h of the shell (see Figure Q}. The stress 
resultants and couples obtained in this way are defined per unit of arc length on the reference surface: 

fh/2 

Ni = / & U (1 + C/R)d( (i = 9,<p), 

J-h/2 

fh/2 

= / &o v (1 + C/R) d(, 

J-h/2 

fh/2 

Qi = / (i + C/R) d( (i = 8,<p), (17) 

J-h/2 
fh/2 

Mi = / & H (l + C/R)CdC (i = 

J-h/2 

t-h/2 



M 6tp = 



M v g = \ & Bv (l + C/R)t<K- 

J-h/2 



h/2 

We evaluate these integrals with the expansion (fT6|) , The tangential stress resultants are 

Ng = K(e g + ve° v ) , 

N v = K(e v + ve°g) , (18) 
Explicit expressions for the transverse shearing stress resultants Qi are not needed since these quantities 
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will be determined from the equilibrium equations. The moment resultants are 

1 



(19) 



2 V # 



7e v 



The extensional rigidity K and the bending rigidity D are defined by 

Eh 



K = 
D = 



Their dimensionless ratio £ is a large number, 



1-J/ 2 ' 

Eh 3 

12(1-J/ 2 ) 



X 12i? 2 



the inverse of which will serve as an expansion parameter for thin shell theory. 



(20) 
(21) 

(22) 



2.2.4 Equilibrium equations 

With the metric JHJ), the components 9, ip and ( of the equilibrium equations {7} respectively become 
(R + C) ((sm8aee) e + &e<p,<p ~ cos6>ct w + sinflo^J + sin (9 [{R + () 2 = 0, 

(R + C) (^{sm9a-0 V ) g +a 9l p >lp + cos9a 6v + sin 03^ J + sin0 (^(R + () 2 a^ v ^j ^ = 0, 

(R + C) ((sin 6>o- ec ) e + a V £, v - sin 8 (& ee + <7 W )J + sin0 {{R + £) 2 &CCJ = °i 



where the equations have been multiplied by sin 9 (R + £), (i? + £) and sin 6* (i? + £) 2 , respectively. The 
stress components are given in a normalized basis. 

The integration on £ of these three equations in the range [—h/2, h/2] yields the equilibrium equations 
for the forces: 

(sin 9 N e ) e + Ne VtV - cos 9N V + sin 9Qg + R sin 9q e = 0, (23) 
(sm6Ng v ) e +N lpilp + cos9N eip + sin9Q v +Rsindq v = 0, (24) 
(sm9Q e ) je + Q lp , lp -$in9 (N e + N tp )-Rsm9q = 0, (25) 
where qg and q v are the components of the tangential load vector per unit area of the reference surface: 



(R + C) a Ci 



h/2 
-h/2 



R 2 q % 



(i = 9,<p) 



We choose the convention that tensile stresses are positive (see Figure [T]). The transverse load per unit 
area of the reference surface is noted q and is taken to be positive toward the center of the sphere: 



h/2 
-h/2 



-R 2 q. 



(26) 



The first two equilibrium equations for the stresses can also be multiplied by C before the integration 
to yield the equilibrium equations for the moments: 



(sinflMfl) e + Mg 9tV - cosflM^ - R sin9Qe = 0, 
(sin 9 M$ v ) e + M v<v + cos 9 Mg^ — R sin 9 Q v = 0. 



(27) 
(28) 



We have neglected small terms in 



h/2 
-h/2 



where i = (9, ip). A third equilibrium equation for 



the moments exists but has the form of an identity: Mg v = M v g. 
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3 Resolution 



3.1 Available methods 



At this stage the elastic theory for a thin spherical shell involves 17 equations: six strain-displacement 
relationships (TilI| - [T4")) . six stress-strain relations (ri"5|) - (|T9l making Hooke's law, and five equilibrium equa- 
tions (f2"3"f - (|23p and (f?7)) - (j2"5)) . There are 17 unknowns: six strain components (eg,e° 



o o o 0\ 

fc y> Wipi "'6! "v> ' ) 



three displacements (w,Vg,v v ), three tangential stress resultants (Ng, N v , Ng v ), two transverse shearing 
stress resultants (Qg,Q v ), and three moment resultants (Mg,M v ,Mg ip ). The three equations (fl"6|) are 
also needed if the tangential stresses (crgg,crg v ,a vl p) are required. The quantities of primary interest to 
us are the transverse de flec tion and the tangential s tresse s (sometimes tangential strain is preferred, as 

We thus want to find the minimum set of 



Sandwell et al 



1997| or 



Banerdt and Golombek 



200 



3)- 



equations that must be solved to determine thes e quantities . 



We are aware of two methods of resolution 



Novozhilov 



1964 



p. 66]. In the first one, we insert the 



strain-displacement relationships into Hooke's law, and substitute in turn Hooke's law into the equilibrium 
equations. This method yields three simultaneous differential equations for the displacements. Once the 
displacements are known, it is possible to compute the strains and the stresses. The seco nd method 
supplements the equilibrium equations with the equations of compatibility NovozhilotA , Il964l . p. 27] that 



relate the partia l derivatives of the strain components. It is then convenient to introduce the so-called 
stress functions Krausl , 119671 p. 243], without direct physical interpretation, which serve to define the 



stress resultants without introducing the tangential displacements. Equations relating the transverse 
displacement and the stress functions are then found by applying the third equation of equilibrium and 
the third equation of compatibility (it is also possible to use all three equations of compatibility in order 
to directly solve for the stress and moment resultants). Once the transverse displacement and the stress 
functions are known, stresses can be computed. 

If the shell thickness is constant, the deformations of a thin spherical shell can be completely calcu- 
lated with both methods. If the shell thickness is variable, the three equations governing displacements, 
obtained with the first method, cannot be decoupled and are not easy to solve. Kraus' method with 
stress functions leads to a system of three equations (relating the transverse displacement and the two 
stress functions), in which the first equation is decoupled and solved before the other two. This method 
thus provides a system of equations much easier to solve and will be chosen in this article. 

When solving the equations, one usually assumes from the beginning the large £ limit, i.e. 1 +£ = £ 
where £ is defined by equation (|22|) . We will only take this limit at the end of the resolution. This 
procedure will not complicate the computations, since we have to compute anyway many new terms 
because of the variable shell thickness. 



3.2 Differential operators 

We will repeatedly encounter the operators 0, which intervene in the expressions (Ti~4"]) for the flexural 
strains (Kg, K®, r ). Since we are looking for scalar equations, we need to find out how the operators Oi 
can be combined in order to yield scalar expressions, i.e. expressions that are invariant with respect to 
changes of coordinates on the sphere. 

The first thing is to relate the Oi to tensorial operators. Starting from the covariant derivatives on the 
sphere Vj, we construct the following tensorial differential operators of the second degree in derivatives: 

n,, = ViV,- + 5 y , (29) 

where Vj denotes the covariant derivative (see Appendix 17. 2\ . These operators give zero when applied 
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on spherical harmonics of degree one (considered as scalars) : 



V ij Y lm = (m = -1,0,1). (30) 

This property can be explicitly checked on the spherical harmonics (| 109[) with the metric and the formulas 
for the double covariant derivatives given in Appendix 17.41 

In two-dimensional spherical coordinates (#, tp), the three operators 0j defined by equations (|15p are 
related to the operators T>ij acting on a scalar function / through 

Oif = V ee f, 

2 f = cac 2 9V w f, (31) 
3 / = csc6V ev f. 

The operators 01,2,3/ actually correspond to normalized T>ijf, i.e. X>y / /{^/gugjj)- 

The usual derivatives of the operators Oi satisfy the useful identities (|112[) - (|113[) which are proven 

in Appendix 17.71 These identities are the consequence of the path dependence of the parallel transport 

of vectors on the curved surface of the sphere. 

Invariant expressions are built by contracting all indices of the differential operators in their tensorial 

form. The indices can be contracted with the inverse metric or with the antisymmetric tensor e 13 (see 

Appendix 17. 4p . which should not be confused with the strain tensor e^. In the following, a and b axe 

scalar functions on the sphere. 

At degree 2, the only non-zero contraction of the T>ij is related to the Laplacian (|104D : 

A'a = g %3 T>ij a 

= (A + 2) a (32) 
= (0i+0 2 )a 

= a,e,e + cot 9 a,e + esc 2 9 a, v , v + 2 a . 

At degree 4, a scalar expression symmetric in (a, 6) is given by 

A(a ; b) = [A'a] [A'b] - [V tJ a] [V 13 b] , 

= [A a] [A 6] - [ViV^a^VH] + [Aa]b + a[Ab] +2ab (33) 
= [Ox a] [0 2 b] + [0 2 a] [0! b] - 2 [0 3 a] [0 3 b] 

= (a,e,0 + a) (esc 2 9 b^^ + cot 9 b t $ + b) + (esc 2 9 a tVtV + cot 9 a <? + a) (b t e,e + b) 
—2 esc 2 9 (a — cot 9 a y ) (b t $ >ip — cot 9 b >v ) . 

where upper indices are raised with the inverse metric: V 13 = g lk g 3l T>ki- The action of an operator does 
not extend beyond the brackets enclosing it. 

If a is constant, A{a ;b) — a A'b. It is useful to define an associated operator Aq that gives zero if 
its first argument is constant: 

Ao(a;b) = A(a;b) - a[A'b] . (34) 
A scalar expression of degree 4 antisymmetric in (a, 6) is given by 

B x (a; b) = g 13 e kl \V lk a] [V jt b] 

= g 13 e kl [V l V k a}[V J Vib] (35) 

= [(0! - 2 ) a] [0 3 6] - [0 3 a] [(0! - 2 ) b] 

= esc 9 {a,fi,8 — esc 2 9 a, v , v — cot 9 a y g) (b,e,ip — cot 9 b jV ) 



CSC 



9 (fl,e,ip — cot 9 a tV ) {b t e,B — esc 2 9 b tVtV — cot 9 6g 
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We will also need another operator of degree 4: 

B 2 {a;b) = [V i( z HVjA'fc] 

= csce(a, e [b!b\ v -a v [b!b}^ . (36) 

The sum of the operators B\ and B 2 is noted B: 

B(a;b) = B 1 (a;b) + B 2 (a;b). (37) 

If either a or b is constant, B(a ; b) = 0. 

The operators A and B have an interesting property proven in Appendix 17.81 for arbitrary scalar 
functions a and b, A(a ; b) and B{a ; &) do not have a degree-one term in their spherical harmonic expansion. 
This is not true of Aq, B\ and B 2 . 



3.3 Transverse displacement 

3.3.1 Resolution of the equations of equilibrium 



The first step consists in finding expressions for the moment resultants (Mg, M v , Mg v ) in terms of the 
transverse displacement and the stress resultants. The extensional strains (e$, e° ,7^) can be eliminated 
from the equations for stress and moment resultants (|18[) - (fTO)) . The flexural strains (k@, k°,t°) depend 
on the transverse displacement w through equations (fT4)) . We thus obtain 



M v = -^(0 2 + u0 1 )w + ^N v , (38) 



D R 
M 9v = -^(1-v)0 3 w + jN 6 



where the parameter £ is defined by equation (|22|) . 

The second step consists in solving the equilibrium equations for moments in order to find the 
transverse shearing stress resultants (Qg,Q v )- We substitute expressions (l38|) into equations ([27]) -([28 |) . 
Knowing that the stress resultants satisfy the equilibrium equations (|23[) - (124[) . we obtain new expressions 
for Qg and Q v (identities (| 1 12|) - f| 1 13[1 are helpful): 

Qe = -^{D\'w) fi -{l-rj)Rq e 

+ 4? (1 - v) (D, e 02 - esc 9 D :ip 3 )w-- {ij.e N g + esc 6 n „ N 9ip ) , (39) 
H 77 

Q v = -■jjgC8cB(DA'w) >ip -(l-ri)Rq v 

+ -£s(l-v) {-D.0 O3 + esc 8D^Oi)w-- (rj.g N 9ip + esc 9 n v N v ) . (40) 

H T] 

The operator A' is defined by equation (f3"2")) and r\ is a parameter close to 1: 

£ ( h 2 Y 1 



i) 



The third step consists in finding expressions for the stress resultants (Ng, N v , N 9lp ) in terms of 
stress functions by solving the first two equilibrium equations (|23|) -(|24 |) . Let us define the following linear 
combinations of the stress and moment resultants: 

(P 8 ,P v ,Pg v ) = (n 6 + ^Mg, N v + ^M v , Ng v + ^M 9 A . (42) 
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We observe that these linear combinations satisfy simplified equations of equilibrium: 

(sin9 Pg) g + Pg v ^ ~ cos6 P v + R sin9 qg = 0, 
(sia9Pe v ) e + P v , v + cos9Pg v + R sin (9^ = 0. (43) 

Comparing these equilibrium equations with the identities (|112p - (|113p . we see that the homogeneous 
equations (i.e. equations (|43[) with a zero tangential load qg = q v = 0) are always satisfied if 

(Pe,P v ,Pe v ) = (0 2 ,0 1 ,-0 3 )F, (44) 

where F is an auxiliary function called stress function. For the moment, this function is completely 
arbitrary apart from being scalar and differentiable. 

Particular solutions of the full equations (|4"3")) can be found if we express the tangential load q T = 
qg9 + q v ip in terms of the surface gradient of a scalar potential £1 (consoidal or poloidal component) and 
the surface curl of a vector potential Vr (toroidal component): 

q T = _Iyr! + |vx (Vr). (45) 

Surface operators are defined in Appendix 17.51 where the terms consoidal/poloidal are also discussed. 
The covariant components of are (qg, sin9q ip ) and can be expressed as — + -^g^^ikVj, which 
gives 

qe = ~li n - 9 + R-ke v ^> 

1 sin 9 

sm9q v = --n, v --—Vg. (46) 
H K 

If the tangential load is consoidal (V = 0), a particular solution of equations (|43|) is given by 

(P e ,P v ,Pe v ) = (1,1, 0)ft. (47) 

If the tangential load is toroidal (Q = 0), a particular solution of equations (|43|) is given by 

(Pg,P v ,Pg v ) = (20 3 ,-20 3 ,0 2 -0 1 )H, (48) 

where we have introduced a second stress function H which satisfies the constraint 

A'H = -V + V Q , (49) 

where Vb is a constant (identities (|112[) - (|113p arc useful). This equation allows us to determine the stress 
function H if the toroidal source V is known. 

The general solution of the equations (|43|) is given by the sum of the general solution ((44)) of the 
homogeneous equations and the two particular solutions (|4"T1) - (|4"8")) of the full equations: 

(Pg, P v , P 9tp ) = (0 2 F + n + 20 3 H, OiF + Sl- 20 3 H, -0 3 F + (0 2 - Ox) H) . (50) 

The stress resultants (Ng, N^, Ng v ) can now be obtained from (Pg, P v , Pg v ) by using equations (|42|) and 

(N g , N v ,N 9v ) = v (P e , P v , Pg v ) + v ^ (O x + v0 2 ,0 2 + vO x ,(\ - v) 3 ) w , (51) 



which finally give 



Xo = 'l[0 2 F + j^ (A' - (1 - v) 2 ) w + fi + 2 OJI j . 
V = ii(OxF + ^(A' -(l-u)0 1 )w + n-20 3 ll ) . (r,2) 
AV = // ( -0 3 F+(l-u)^0 3 w+(0 2 -0 1 )H 
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The fourth step consists in expressing the third equation of equilibrium (j25|) in terms of the transverse 
displacement w and the stress functions F and H. For this purpose, it is handy to express the transverse 
shearing stress resultants (Qe,Qip) in terms of (w,F,H). We thus substitute N v , N$, Ng v , given by 
equations (j52")) . into the expressions for Qg and Q v , given by equations (p?9"|) - (|4T))) : 

Qe = -^{riD6!w) )0 + ^^{^D) fi O 3 -aac9{fiD) )V O^w 

- (V,e 2 - esc 9 r hv O 3 )F+n :0 - (rjtyj - esc 9 V v 

- 2 (r?, e 3 + esc 9 V , v 2 ) H + csc9 (r hv A'H - n{A'H) >v ) , (53) 

Q v = -^csc9 (r)DA'w) tlp + ^^ (- ( V D) e 3 + csc9 ( V D)^0^ w 

- (-V,e0 3 + esc 9r), v Gi) F + esc 0(O i¥ , - (r]Q), v ) + Vg 

+ 2(r lt e0 1 + csc9T], ¥> 3 )H-(r] i eA'H-r ] (A'H),e) . (54) 

The various terms present in and Q v can be classified into generic types according to their 
differential structure. Each generic type will contribute in a characteristic way to the third equation 
of equilibrium. It is worthwhile to compute the generic contribution of each type before using the full 
expressions of the stress resultants with their multiple terms. Identities (|112[) - (|113[) arc helpful in this 
calculation. We thus write the third equation of equilibrium (f2"5"j) as follows: 

l(Qe;Q v ;0)-(N e + N v )-Rq = 0, (55) 

where the differential operator X is defined for arbitrary expressions (X, Y, Z) by 

1(X;Y;Z) = esc 9 (sin 6 X) g + esc 9Y V - cot 9 Z. (56) 

This operator must be evaluated for the following generic types present in (Qg, Q v ): 











I{a,e 


esc 9 


0) 


= Aa, 










I (— CSC 


9 a i¥ , ; a t e 


0) 


= o, 




o 2 b 


— csc 9 


3 b;~a fi 


3 b + esc 


9a, v 1 b 


0) 


= A (a;b), 




o 3 b 


+- esc 9 


2 b;-a fi 


Oib~ esc 


9a, v 3 b 


0) 


= B x (a; b), 


1( 


- CSC 


9(a v A'b 


-a[A'b] tV 


);a e A'b- 


a[A'b} t e 


0) 


= 2B 2 (a;b), 



where a and b are scalar functions. The operators Ao, B\ and B 2 are defined in section f3. 21 

We now substitute (Qe, Qip) and (N v , Ng) into the third equation of equilibrium ([55]) and use formulas 
(f57|) . We thus obtain the first of the differential equations that relate w and the stress functions F and 
H: 

A' (r)D A'w) ~{l-v) A(r)D ;w) + R 3 A(r) ; F) + 2R 3 B{q ; H) = -i? 4 q + R 3 (AQ - A'(r)Cl)) . (58) 
The operators A', A and B are defined in section [5721 

3.3.2 Compatibility relation 

A second equation relating the transverse displacement and the stress functions comes from the compat- 
ibility relation which is derived by eliminating (vg,v v ) from the strain-displacement equations (|13[) : 

(sin 9 7 °^) g = (sin 2 9 e%) g + e° e<tp<tp - sin 9 cos 9 e% + 2 sin 2 9 e° e - ^ A'w . (59) 

The strain components are related to the stress resultants through Hookc's law (fT8|) . The substitution of 
equations (fT5|) into the compatibility equation (fBT))) gives 

A'(a(N e + N ip ))--A'w-(l + v)J(aN e ;aN ip ;aN e(p ) =0. (60) 
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where a is the reciprocal of the extensional rigidity: 

1 1 



(61) 



K{1 - v 2 ) Eh ' 
For arbitrary expressions (X, Y, Z), J(X ;Y ;Z) is defined by 

J(X;Y;Z) = esc 2 9 ((sin 2 9X, ) fi +Y, V)V + 2 (sin (9 Z v ) e ) - cot 8 Y fi + 2 Y . (62) 

As in the case of the third equation of equilibrium, it is practical to classify the terms present in 
(N$, N v , Ngtp) into generic types and evaluate separately their contribution to the equation of compati- 
bility. There are three types of terms for which we must evaluate the operator J (identities (| 11 2|) - (| 11 5[) 
are helpful in this calculation): 

J(a;a;0) = A' a, 
J(a0 2 b;a0 1 b;-a0 3 b) = A(a;b), (63) 
J{2a0 3 b;-2a0 3 b;a{0 2 -0 1 )b) = 2B(a;b), 



where a and b are scalar functions. The operators A', A and B are defined in section [3~2 
We now evaluate J in equation (|60p with expressions (|52[) and formulas 



J(aNe;aN 9 ;aNe v ) = A(rja ; F) + A' (rjaD A'w) - (1 - v) A{r)aD ;w) 

+ A'(j]aCL) +2B(r]a;H). (64) 
The term A(r)aD ; w) can be rewritten with the help of the following equality: 

1 - v 2 

„ AOqaD ; w) = A w - A(n ; w) , (65) 

since (1 — v 2 )r\aDjR 2 — rj/£ and (?7/0,i = — V,i- 

We finally substitute A^g + N v , given by equations (|52|) . into the compatibility equation (|60|) and 
use expressions ([64| -([65 |l . We thus obtain the second of the differential equations that relate w and the 
stress functions F and H: 

A' (rja A'F) -(l + v) A{r)a ;F) — — A(rj ; w) - 2(1 + v) B(r]a ;H) = -(1 - v) A'(rja 0) . (66) 

R 

3.4 Tangential displacements 



Assuming that the flexure equations (f58|) and (|66|) for the transverse displacement w and the stress 
functions {F, H) have been solved, we now show how to calculate the tangential displacements. 

In analogy with the decomposition of the tangential load in equations (|4"5"1) - (|4"6")) . the tangential 
displacement can be separated into consoidal and toroidal components: 



v = VS + Vx(Tf), (67) 

where S and T are the consoidal and toroidal scalars, respectively. The covariant components of v are 
(ve, sin 9 v v ) and can be expressed as + gi k £ikT,j (see Appendix [73]), which gives 

v g = S t e + esc 6 T j(p , 
sm6v v = S^-wa.BT fi . (68) 
The strain-displacement equations (|13[) become 

4 = ^((0 1 -l)S + 3 T + w) , 

e° = ^((0 2 -l)S-0 3 T + w) , (69) 
lg v - ^(2 3 S+(0 2 -0 1 )T) . 
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The stress resultants (Ng, N v , Ng v ) given by equations (1181) become 

N e = ^-(AS+(l + is)w-(l-v)((0 2 -l)S-0 3 T)) 
K 



K 
~R 

K 1 - v 



N v = — (A5+(l + i/)tfl-(l-i/)((Oi-l)5 + 3 T)) , (70) 



N 0V = ___(20 3 5+(02-Oi)r) . 
The toroidal potential T cancels in the sum iVg + N v : 

N e + N v = ^ (1 + !/) (A S + 2tu) . 
rt 

The consoidal displacement potential 5 can thus be related to (w, F, f2) by using expressions (|52"|) for the 
stress resultants: 

AS = Ri]a (1 - v) (A'F + 2 fi) + | A'to - 2w , (71) 

where a is defined by equation (|STjl . 

It is more difhcult to extract the toroidal displacement potential T. When the shell thickness is 
constant, decoupled equations for the displacements can be found by suitable differentiation and combi- 
nation of the three equilibrium equations (|23|) - ((2"5|) for the stress resultants. This method does not work 
if the shell thickness is variable because the resulting equations are coupled. The trick consists in relating 
the tangential displacements to (w, F, H) by the way of equations similar to the homogeneous part of the 
first two equilibrium equations, but with (Ng, N v , Ng v ) replaced by ^(Ng,N v , Ng v ), so that derivatives 
of j£ do not mix with derivatives of T. We will thus calculate the following expression in two different 
ways (from equations (f52"|) and ([70]) ): 



Z ( -| N e ; -| N v ; ^ Ng v ) = esc 2 t) ( X. , sin H V,,) . 



where X and Y are defined by 

X = sinOX (^Ng;^N ev ;^N lp i . 

Y = S mei(^N eip ;^N v ;~N 6 

The operator I is defined by equation (|5S|) . 

As before, it is easier to begin with the evaluation of the operator Z for generic contributions: 

Z(a;a;0) = 0, 
Z(a0 2 b;a0 1 b;-a0 3 b) = -B(a;b) , 
Z(2aO s b ; -2a(D 3 b ; a(0 2 - Oi)b) = 2A(a;b)- A' (aA'b) , 

On the one hand, the evaluation Z for (Ng,N v , Ng v ) given by equations ([70|) gives 

1 V AA'T . 



2 

On the other hand, the evaluation of Z for (Ng, N v , Ng^) given by equations ([52")) gives 

-Rb(^;F)+(1-u)B^; W S )+2Ra(^;H)-RA'(^A'h) . 

The equality of the two previous formulas yields the sought equation for T: 

AA'T = 2R (1 + v) (B (r)a ;F) — 2A (rja ;H) + A' (rja A'H)) + 2 B (r] ; w) . (72) 
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We have used the relation (?7/£),i = — V,i- 

The equation for T shows that toroidal displacement always occurs when the shell thickness varies. 
The right-hand side of equation (172]) only vanishes when two conditions are met: (1) there is no toroidal 
source (so that the stress function H vanishes) and (2) the shell thickness is constant (so that the terms 
in B vanish). 

3.5 Stresses 

Stresses can be computed from (w,F,H) and f2 by substituting equations (fT4")h (fT5)) and ([52"]) into equa- 
tions (HU): 



n v , o o * tt\ , E fv C 

R(l-v 2 ) V? R + C 



x w = -L(0 1 F + n-20 3 H) + — —[j~^—:)(A'~(l~iy)0 1 )w, (73) 



Stresses at the surface are obtained by setting £ = /i/2. 

4 Flexure equations and their properties 
4.1 Thin shell approximation 

The flexure equations derived in section [3] already include several assumptions of the thin shell theory, 
but not yet the first one stating that the shel l is thin. Of course, the three other assumptions can be 



seen to be consequences of the first one [see \Kraus\ . I1967L chap. 2.2], but we have not imposed in a 



quantitative way the thinness condition on the equations. We thus impose the limit of small h/R or, 
equivalently, the limit of large £ (defined by equation (|2"2"f ) on the flexure equations for (H, F,w, S,T), 
This procedure amounts to expand r) « 1 — l/£ (neglecting terms in l/£ wherever appropriate) and to 
neglect the derivatives of ij in equations (|4U)) . (|58p . (|6"6")l . (fTTj) and ([72")) . We thus obtain the final flexure 
equations for the displacements (w,S,T) and for the stress functions (F,H): 

A'H = -V + V a , (74) 

E>3 

A' (D A'w) -(l-v) A(D ; w) + R 3 A'F = -R 4 q-2 R 3 (l + — Ail , (75) 

A' {a A'F) -(1 + v) A(a ; F) — -i A'w; = -(1 - v) A' (a O) + 2(1 + v) B(a ; H) , (76) 

R 

AS = Ra{l-v){A'F + 2tt) + ^A'w-2w , (77) 

A A' T = 2R(l + v) (B (a;F) - 2 A(a;H) + A' (a A 1 ' H)) . (78) 

Recall that the differential operators A', A and B are defined by equations (|3"2"|) . and ([57)1 . The 
potential pairs (f2, V) and (S,T) are related to the tangential load and displacement by equations ([4"5")l 
and (|67|) . respectively. In the second equation, the term i? 3 AJ7/£ has been kept since it could be large if 
Q has a short wavelength. For the same reason, the term A'w/£ has been kept in the fourth equation. 

In the third and fifth equations, the terms depending on H belong to the right-hand sides since they 
can be considered as a source once H has been calculated from the first equation. The same can be said 
of the terms depending on w and F in the fourth and fifth equations: they are supposed to be known from 
the simultaneous resolution of the second and third equations. The difficulty in solving the equations 
thus lies with the two flexure equations for w and F which are linear with non-constant coefficients (all 
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other equations are linear - in their unknowns - with constant coefficients) . Once these two core equations 
have been solved, all other quantities are easily derived from them. 

The bending rigidity D, defined by equation (|21[) . characterizes the bending regime: the shell locally 
bends in a similar way as a flat plate undergoing small deflections with negligible stretching. The pa- 
rameter a, defined by equation (|6ip . is the reciprocal of the extensional rigidity K and characterizes the 
membrane regime of the shell in which bending moments are negligible and the load is mainly supported 
by internal stresses tangent to the shell. Since D and K are respectively proportional to the third and 
first power of the shell thickness, the membrane regime (in which the D-depending terms are neglected) 
is obtained in the limit of an extremely thin shell. This observation and the fact that such a shell, lacking 
rigidity, cannot support bending moments justify the use of the term 'membrane'. 

The weight of the various terms in the equations depends on two competing factors: the magnitude 
of the coefficient multiplying the derivative and the number of derivatives. On the one hand, a coefficient 
containing D will be smaller than a coefficient containing oT x since aD/R 2 ~ is a small number 
(see equation ). On the other hand, a large number of derivatives will increase the weight of the term 
if the derived function has a small wavelength. The transition between the membrane and the bending 
regimes thus depends on the wavelength of the load: if the load has a large wavelength (with respect to 
the shell radius), the flexure of the shell will be well described by equations without the terms depending 
on D (see section UTTj) . whereas the flexure under loads of small wavelength is well described by equations 
keeping only the D-depending terms with the largest number of derivatives (see section f5.2[) . 

Stress functions are associated with membrane stretching and give negligible contributions in the 
bending regime. Formulas (|52[) show that F (respectively H) plays the role of potential for the stress 
resultants in the membrane regime when the load is transversal (respectively tangential toroidal) . There 
is no stress function associated with the tangential consoidal component of the load because it is identical 
to O, the consoidal potential of the load. 

The variation of the shell thickness has two effects. First, it couples the spherical harmonic modes 
that are solutions to the equations for a shell of constant thickness. Second, the toroidal part remains 
intertwined with the transversal and consoidal parts, whereas it decouples if the thickness is constant. 
For example, the toroidal load is a source for the transverse deflection through the term B in the third 
flexure equation. Furthermore, the stress function F is a source for the toroidal displacement in the fifth 
flexure equation. 

For numerical computations, the flexure equations (|74j) -(|78 |l obtained in the thin shell approximation 
(to which we can add the equations (|73|) for the stresses) are adequate. For this purpose, it is not useful 
to keep small terms in l/£ since the theory rests on assumptions only true for a thin shell. In the rest of 
the article, we will continue to work with the equations (|49]) . (|58|). (|66|) . (|7T|) and (|72)l for more generality. 

4.2 Covariance 

Because of their tensorial form, the flexure equations ([74|) -(f78 f are covariant; this is also true of the more 
general equations P5|) , ([55]) . (fTTj) and (fT2")l . This property means that their form is valid in all 

systems of coordinates on the sphere, though the tensor components (the covariant derivatives of scalar 
functions, the metric and the antisymmetric tensor) will have a different expression in each system. The 
scalar functions will also have a different dependence on the coordinates in each system. The covariance 
of the final equations was expected. We indeed started with tensorial equations in section 12. l\ their 
restriction to the sphere in principle respected the tensoriality with respect to changes of coordinates 
on the sphere. However the covariance of the two-dimensional theory was not made explicit until we 
obtained the final equations. This property is thus a strong constraint on the form of the solution and 
a check of its validity, though only necessary and not sufficient (other covariant terms may have been 
ignored) . 
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Another advantage of the covariant form is the facility to express the final equations in different 
systems of coordinates (even non-orthogonal ones) with the aim of solving them. For example, the fi- 
nite difference method in spherical coordinates (#, <p) suffers from a very irregular grid and from pole 
singu larities. These problems can be avoided with the 'cubed sphere' coordinate system 



Ronchi et al 



1996] . Operators including covariant derivatives (here the Laplacian and the operators A and B) can ex- 



pressed in any system of coordinates whose metric is known. Christoffel symbols and tensorial differential 
operators can be computed with symbolic mathematical software. 



4.3 Degree one 

Displacements of degree one require special consideration. All differential operators acting on (w,F,H) 
(that is A', A and B) in the flexure equations (f58"| and (|66p can be expressed in terms of the operators 
T>ij = V^Vj + gij (see section [372)1 . The Z>y have the interesting property that they give zero when acting 
on spherical harmonics of degree one (see equation ([30]) ). Therefore the degree-one terms in the spherical 
harmonic expansion of w vanish from the flexure equations. The magnitude of the transverse deflection 
of degree one neither depends on the load nor on the elastic properties of the spherical shell. 

More generally, the homogeneous (q = Q = V = 0) flexure equations (PE))) . (155)1 and (|6"6"j) are satisfied 
by w and F being both of degree one. According to equations (f7T)) - (f72|) . the corresponding tangential 
displacement is constrained by AS = —2w and AA'T = 0, so that S and T are also of degree one, 
with S = w. These conditions lead to vanishing strains (see equations (|69p ) which indicate a rigid 
displacement. The total displacement is then given by 

u = w r + Vw + V x (Tf ) , 

where w and T are of degree one. In Appendix 17. 6[ we show that the first two terms represent a 
rigid translation whereas the last term represents a rigid rotation. As expected, stresses vanish for 
such displacements (see equation ([73]) ). This freedom in translating or rotating the solution reflects the 
freedom in the choice of the reference frame (in practice the reference frame is centered at the center of 
the undeformed shell). The same freedom of trans lation is also found in the theory of de f ormations of a 



1997 



Blewitt 



spher ical, radially stratified, gravitating solid [e.g. \Farrell Il972t I Greff-Lefftz and Learos 

hooaj . 

What can we say about degree-one loading? Let us first examine what happens when the flexure 
equations are projected on the spherical harmonics of degree one. Since the operator A' annihilates the 
degree one in any spherical harmonic expansion, terms of the form A'/ vanish when they are projected 
on the spherical harmonics of degree one. Moreover, the operators A and B also vanish in this projection 
since they do not contain a degree-one term (see equations fjl 16|l - f)117[) of Appendix 17. 8p . Therefore the 
degree-one component of the flexure equation (|66p is identically zero, whereas the degree-one component 
of the flexure equation (|58p is 

/ dw(? + V.q T )y i = (i = x,y,z), (79) 
Js 

where cLo = sin 9 d9 dip and Y; are the real spherical harmonics of degree one. The integral is taken over 
the whole spherical surface. We have used the relation Af2 = — R V • Qt derived from equations ([35)1 and 

(fim 

The first term in the integrand of equation (|T9"|) is the projection on the Cartesian axes (x, y, z) of 
the vector field q f : 

(qY x ,qY y ,qY z ) = (qr ■ x, qr ■ y, qr ■ z) , 
where we have used formulas (|109p for the spherical harmonics. 
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The second term in equation (|79|) can be rewritten with the identity () 1 0211 and Gauss' theorem (|107j) : 

/ dw (V • q r ) Yt = — du> qx • V Yi (i = x, y, z) , 
Js Js 

where the Yi are considered as scalars. Since qr is orthogonal to f , the integrand is the projection of qr 
on the Cartesian axes (x,y,z): 

(q T • V Y x , q T • V Y y , q T ■ V Y z ) = (q T ■ x, q T ■ y, qr ■ z) , 

where we have used formulas (111 0[) for the gradients of the spherical harmonics. 

Recalling that the transverse load q was defined positive towards the center of the sphere (see equation 
(l2l)|) h we define a total load vector q = — qr + qr- With the above results, we can rewrite the degree-one 
projection (|79j) as 

/ dw(q-x,q-y,q- Z ) = (0,0,0), (80) 
Js 

which means that the integral (over the whole spherical surface) of the projection on the coordinate axes 
of the total load vector vanishes. This result is the consequence of the static assumption in the equations 
of motion fT|). since a non-zero sum of the external forces would accelerate the sphere. 



In practice, degree-one loads on planetary surfaces are essentially due to mass redistribution Greff-Lefftz and Leans 



19971 ] and have a tangential consoidal component (for example the gravitational force is not directed to- 
ward the center of figure of the shell). If the shell thickness is variable, a non-zero f2 of degree one will 
induce degrees higher than one in w and in S. If the shell thickness is constant, the degree-one load 
drops from the flexure equations and w is not affected. However, the degree-one f2 generates 

(assuming a constant shell thickness) a degree-one tangential displacement through equation (|7ip , so that 
S w. Whether the shell thickness is variable or not, a degree-one 17 thus genera tes a total displacement 



which is not only a translation but also a tangential deformation 
do not vanish as shown by equations (I73p . 



Blewitt 



in which case stresses 



5 Limit cases 
5.1 Membrane limit 

A shell is in a membrane state of stress if bending moments {Mg,M Vl M@ v ) can be neglected, in analogy 
with a membrane which cannot support bending moments. Equations (|19p show that this is true if the 
bending rigidity vanishes: D — 0. Consistency with equation (|2"2"|) imposes the limit of infinite £ (or 
T) = 1). With these approximations, the first flexure equation ([58"]) for (w,F,H) becomes 

A'F = -Rq-2Q. (81) 

The second flexure equation for (w, F, H) becomes 

- A'w = A' (a A'F) - (1 + v) A{a ; F) - 2(1 + v) B(a ;H) + (l-u) A' (a ft) . (82) 
it 

If q is independent of w, F and w can be successively determined with spherical harmonic transforms 
from equations (|8Tj) and ([82]) (though the right-hand side of the latter equation must be computed with 
another method). If q has a linear dependence in w (such as when the sphere is filled with a fluid), w 
can be eliminated between equations (f8~T|) and (l82|) . so that F and w can also be computed in succession 
(however the equation for F cannot be solved by a spherical harmonic transform). H is supposed to be 
known since equation (|19")) is not modified and can be solved with spherical harmonics. 
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The stresses are obtained from equations ([73]) with the additional approximation of neglecting the 
term in £: 

1 

h 



{aee^w, &e v ) = r (°i F + Q + 2 °3#, O x F + fi - 20 3 H, -O a F + (0 2 -Ox)H) 



Bending moments play a small role if the load has a large wavelength. In practice, the threshold at 
which bending moments become significant can be evaluated from the constant thickness equation for w 
(see section f5 . 3 . 1 [) . One must compare the magnitudes of the terms in D and I /a in the left-hand side 
of equation (|88|) . Bending moments are negligible (i.e. the term in D) if the spherical harmonic degree I 
of the transverse displacement w is such that 

*s fc yf> (83) 

where k = (12(1 - v 2 )) 1 ' 4 w 1.8 if we take v = 1/4. This threshold is about 10 for a planet with a radius 
of 3400 km and a lithosphcric thickness equal t o 100 km. Flexure equat ions in the membrane limit of a 
shell with constant thickness have been used by Sleev and Phillivs 1985| to study the lithospheric stress 
in the Tharsis region of the planet Mars. 

5.2 Euclidean limit 

The equation for the deflection of a recta ngular plate with variable thickness has been derived by 



Timoshenko and Woinowskv-Krieaen 19641 p. 173]. We will check here that the Euclidean limit of 
our equations gives the same answer. 

Let us define the coordinates (x, y) = (R<p, R0'), where 6' = 8/2 — 9 is the latitude. We work in a 
small latitude band around the equator (so that 9'<1) and in the limit of large spherical radius R and 
large £ (rj = 1). Under this change of coordinate, each derivative introduces a factor R so that terms with 
the largest number of derivatives dominate. In particular, covariant derivatives can be approximated by 
usual derivatives. The surface Laplacian (|104p can be approximated as follows: 

A R 2 (— — 

\dx 2 dy 2 

= R 2 A e . 

We assume that there are no tangential loads (f2 = V = H = 0). The flexure equations (|75|) - (f75]) for the 
transverse displacement become: 

R 4 A e (D A e w) — — R 4 A e (D ; w) + R 5 A e F = -R 4 q, (84) 
R 4 A e (aA e F)-{l + v)R 4 A e (a;F)~RA e w = 0, (85) 

where the operator A e is defined by 

A e (a;b) = [A, </• ; A, />) h ,,. 



d 2 a\ fd 2 b\ n( da \( 9 b \ ( d a \ ( 9 h \ 
dx 2 J \ dy 2 J \ dxdy J \ dxdy J \ dy 2 J \ dx 2 J 



Equation (|85[) gives a relation between the magnitudes of F and w: 

0(R 3 A e F) ~ 0(w/a) . 
In the large R limit, equation (|84)) thus becomes 

A e (D A e w) ~{l-v) A e {D;w) = 
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Stark et al 



200 



3], 



Kirbv and Swat 



n 



2004j and 



Perez- Gussinve et al 



which i s the equation der i ved bv\Timoshenko and Woinowskv-Krieaen 11964 1. Th is equation has been 
used by 



2004j for the local analysis 



of the lithosphe re of the Earth. Its one-dimensional version, in which A e vanishes , has been used by 



Sandwell 



1984j to describe the flexure of the oceanic lithosphere on Earth and by 



Stewart and Watts 



19971 ] to model the flexure at mountain ranges. 



5.3 Shell with constant thickness 
5.3.1 Displacements 

If the thickness of the shell is constant, the toroidal part of the tangential displacement decouples. The 
terms in B indeed drop from the flexure equations (|58[) and (|66p so that the equations for (w, F) depend 
only on (q, CI): 



rjD A'A'w -(1-v)t}D A'w + rjR 3 A'F 

A'A'F - (1 + u) A'F - — A'w 
Ra 



R 4 q + R 3 ((1 - rf) Afi - 2rjn) , 

A'n, 



(86) 
(87) 



where we used the property A(a ; b) — aA'b valid for constant a. We eliminate F from these equations 
and obtain a sixth order equation relating w to (q, fl): 

rjD A A'A'w + — A'w = -R 4 (A'-l-v)q+ R 3 ( — !— A' - 1 - v) Afl . (88) 
a ' \ 1 + 6 / 

The elimination of A'F between equations {77]) and (|8^)) gives an equation relating the consoidal 
displacement potential S to (w, q, fl): 



AS 



1 



' AA'w-2w-{l-v)R 2 aq+^—^RaAtl 



1+f 1+C 

Terms not including a Laplacian can be eliminated with equation 
solution for S in terms of (w, q, ft): 

1 



S 



1 + f 1 - v 



' (A + 1 + v) A'w + w + R 2 a q - (A - £(1 +v))Sl, 



1+^ 

so that we obtain an explicit 



(89) 



where the integration constant has been set to zero. 

Equations (|4"9"|) and l|72p give an equation for the toroidal displacement potential: 

A'T = -2r)Ra (1 + v) V , 



(90) 



where the integration constant has been set to zero. We have assumed that AV^ ^ 0, otherwise we get 
A'T = 0. 

The differential equations given in this section can be solved with spherical harmonics so that the co- 
efficients of the spherical harmonic ex pansions of (w. S, T) can be express e d in term s of th e corresponding 
coefficients of the loads (q, il, V) (see 



Kraus 



1967 | 



Turcotte et al 



1981] , 



Banerdt 



19861] ') 



5.3.2 Comparison with the literature 

We now co mpare our equa tions for a shell of const ant thickness w ith those found in the literature. The 
formulas of 



Banerdt 



19861 ] (taken from the work of 



Vlasov 



1964J) are the most general: 



r>2 /i 

£> (A 3 + 4A 2 ) u; + — (A + 2) w = -R 4 (A + 1 - v) q + R 3 ( - A - 1 - v ) AVl , (') l. i 



(A + 2) x = 



R 2 



AV. 



(92) 
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Banerdt's notation is slightly different: his formulas are obtained with the substitutions £ — > ip, fl — ► 
and V — > i2V. The normal rotation \ is proportional to the radial component (in a normalized basis) of 
the curl of the tangential displacement: 

X=^(Vxv),. 
The curl V x v is related to our surface curl (| 1 03[) by 

Vxv = Vxv + esc 6 ((sin9v v ) g — ve,tp^j r . 
With the formulas and (|104p . we get Vxv = — ATfso that equation (|92[) becomes 

A'T = -2Ra{l + v)V . (93) 



We see that Banerdt's equations (|9lj) and (|93|) coincide with our equations l[88]) and (|90|) in the limit 
of large £ (77 = 1), with one exception: the bending term for w is written D(A 3 + 4A 2 )w instead of 
DA A'A'w — (A 3 + 4A 2 + 4A)w. This error has propagated in many articles and is of consequence for 
the degree-one harmonic component, since it violates the static assumption and spoils the translation 
invariance discussed in section |4~51 The impact on higher degrees is negligible. Because of this mistake, 
many authors give a separate treatment to the first harmonic degree. Banerdt also gives formulas for the 
tangential displacements in terms of consoidal and toroidal scalars (A, B) corresponding to our scalars 
(S,T): his formula (A10) is equivalent to our equation (|89|) in the limit of large £. 

If we ignore temperature effects, Kraus' first equation for (w, F) is equivalent to our equation 
(|86|) in the limit of large £, whereas his second equation for (w, F) is equivalent to the combination 



eq.® 



definition of Kraus' stress function F 



Kraus 



eq. (|86p in the limit of large £ [see Krausl 119671 eq. 6.54h and 6.55d]. Note that the 



19671 p. 243] differs from ours: 



^Kraus 



F-k{l-v)^w, 



with k = 1. This freedom of redefining F for arbitrary k remains as long as D is constant. The flexure 
equation for w, equation (|88p . is unaffected so that the solution for w is unchanged. In the final step, 
Kraus makes a mistake when combining the two equations for (w, F) and thus obtains a flexure equation 
for w with the same error as in equation (|9ip . Krau s does not include toroidal loading. The flexure 
equation of 



Turcotte et at 



198l| is taken from 



Kraus 



1967j | without the tangential loading and is the 



same as equation (I^TI) with 
The flexure equation of 



0. 



Brotchie and Silvester 



1969] is given in our notation by 



DA 2 i 



R2 D 4 
— w = —H 

a 



(94) 



where q includes their term jw describing the response of the enclosed liquid. This equation can be 
obtained from our equation (|88|) as follows: keep only the derivatives of the highest order in each term, 
set f2 = 0, take the limit of large £ (rj = 1) and integrate. Brotchie and Silvester choose to work in 
the approximation of a shallow shell and with axisymmetrical loading, solving their equation in polar 
coordinates with Bessel-Kelvin functions. The reduction to fourth order in equation the shallow 

shell approximation and the axisymmetrical assumption are not justified nowadays since the full equation 
(|88p can be quickly solved with computer-generated spherical harmonics. 

The contraction due to a transverse load of degree 0, w = — R 2 a(l — v)q/ 2, is equivalent to the 
radial displacement computed by Love in the limit of a thin shell Lovel 1 1944 p. 142]. However ad- 



ditional assumptions about the initial state of stress and the internal density changes are necessary 
Willemann and Turcotte , 1982| so that the degree is usually excluded from the analysis. 
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5.3.3 Breakdown of the third assumption of thin shell theory 



With the spherical harmonic solutions of the equations for a shell of constant thickness, it is possible to 
check the thin shell assumption stating that the transverse normal stress is negligible with respect to the 
tangential normal stress. The magnitude of the former can be estimated by the load q (see definition 
(|2l)|) ) whereas the magnitude of the latter can be approximated with formulas ([73]) evaluated on the outer 
surface: 



1 . 

2 \ a se + <r vv ) 



— A'F 
2h 



Eh 



4R 2 (1 - v) 



A'i 



where we have assumed the absence of tangential loads (51 = 0) and the limit of large £. We can relate 
(Tt to q by using the solution in spherical harmonics of equations (|8T[) and ([88]). Since the thin shell 
assumption is expected to fail for a load of sufficiently small wavelength, we assume that the spherical 
harmonic degree I is large. Assuming I 3> 1, we obtain 



1 



(A'F) im « — 



Ra 



Wir, 



R 2 a 



1 



i 4 



■ qtn 



where the spherical harmonic coefficients are indexed by their degree £ and their order m. If the shell is 
not in a membrane state of stress (see equation ([531 ). £ 2 > 2R/h so that ot can be approximated by 



Or) 



lm 



4£ 2 



The thin shell assumption holds if q < <7t, that is if 

e< yf^+v)- or 



■ qir, 



A > 



2tt 



(95) 



where A is the load wavelength (A « 2nR/£). We have ^3(1 + u) w 1.9 and 2tt/- x /3(1 + u) w 3.2 if we 
take v — 1/4. This co ndition on A is con s istent with the transition zone betw een the thin and thick shell 



Janes and Melosh 



19901 and 



Willemann and Turcotte 



Zhona and Zuber 



200ot . but does not coincide with 



1982j | . which is I < 2%-^/ R/h (this last condition looks 



responses analyzed in 
the constraint given in 
more like the threshold (f83|) for the membrane regime). 

Though the stress distribution is affected, the limit (|9"5"|) on the degree I is not important for the 
displacements, since they tend to zero at small wavelengths. Therefore the theory does not break down 
at short wavelength if one is interested in the computation of the gravity field associated to the transverse 
deflection of the lithosphere. 



6 Conclusion 

The principal results of this article are the five flexure equations (|74|) - ((78|) governing the three displace- 
ments of the thin spherical shell and the two auxiliary stress functions. Stresses are derived quantities 
which can be obtained from equations J73|) . The shell thickness and Young's modulus can vary, but 
Poisson's ratio must be constant. The loads acting on the shell can be of any type since we extend the 
method of stress functions to include not only transverse and consoidal tangential loads, but also toroidal 
tangential loads. The flexure equations can be solved one after the other, except the two equations 
([75|) -([75 ]) for the transverse deflection w and the stress function F, which must be simultaneously solved. 
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Tangential loading is usually neglected when solving for the deflection because of its small effect. In that 
case, it is sufficient to solve the two equations (| 75 1) - ([76")) with f2 = H = 0: 

A' (D A'w) - (1 - v) A(D ; w) + R 3 A'F = -R 4 q, 

A' (a A'F) - (1 + u) A(a ;F) - — A'w = 0. 

R 



However tangential loading must be taken into account when computing stress fields -B«nenM ll986|. 

In the long- wavelength limit (i.e. membrane regime), all equations can be solved one after the other 
because it is possible to solve for F before solving for the transverse deflection. If a small part of the shell 
is considered, the flexure equations reduce to the equations governing the deflection of a flat plate with 
variable thickness. If the shell thickness is constant, the flexure equations reduce to equations available 
in the literature which can be completely solved with spherical harmonics. Our rigorous treatment 
of the thin shell approximation has clarified the effect of the shell thickness on the flexure equations. 
We emphasize the need to use the correct form for the equations (without the common mistake in the 
differential operator acting on w) in order to have the correct properties for the dcgrcc-onc deflection and 
degree-one load. 

We have also obtained two general properties of the flexure equations. First we have shown that 
there is always a toroidal component in the tangential displacement if the shell thickness is variable. 
Second we have proven that the degree-one harmonic components of the transverse deflection and of the 
toroidal component of the tangential displacement do not depend on the elastic properties of the shell. 
This property reflects the freedom under translations and rotations of the reference frame. Besides we 
have shown that degree-one loads are constrained by the static assumption but can deform the shell and 
generate stresses. 

This article was dedicated to the theoretical treatment of the flexure of a thin elastic shell with 
variable thickness. While the special case of constant thickness admits an analytical solution in terms of 
spherical harmonics, the general flexure equations must be solved with numerical methods such as finite 
differences, finite elements or pseudospectral methods. In a forthcoming paper, we will give a practical 
method of solution and discuss applications to real cases. 
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7 Appendix 

7.1 Covariant, contravariant and normalized components 

Tensors can be defined by their transformation law under changes of coordinates. The two types of 
tensor components, namely covariant and contravariant components, transform in a reciprocal way under 
changes of coordinates. Tensor components cannot be expressed in a normalized basis: the space must 
have a coordinate vector basis (for contravariant components) and a dual basis (for covariant components) 
which are not normalized. 

The only exception is a flat space with Cartesian coordinates, where covariant, contravariant and 
normalized components are identical. Since the metric is the scalar product of the elements of the 
coordinate vector basis, the covariant components arc related to components defined in a normalized 
basis (written with a hat) by 

Ui — Ui , 
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whereas the relation for contravariant components is 

1 



The normalized Cartesian basis (x, y, z) is related to the normalized basis for spherical coordinates 
(r,6,<p) by 

x = cos 9 cos (p 9 — sin tp <p + sin 9 cos ip f , 

y = cos 9 sin tp 8 + cos ip tp + sin 9 sin ip f , (96) 

z = — sin 9 6 + cos 6* f . 

7.2 Covariant derivatives 

Usual derivatives are indicated by a 'comma': 

_ <M_ 
Vi ' j ~ ih-J • 

Covariant derivatives (defined below) are indicated by a 'bar' or by the operator V*: 

v i\j = V j v i ■ 

The former notation emphasizes the tensorial character of the covariant derivative since the covariant 
derivative adds a covariant index to the vector. The latter notation is more adapted when we are interested 
by the properties of the operator. 

The covariant derivative of a scalar function / is equal to the usual derivative, fu = and is itself 
a covariant vector: fu — Uj. Covariant derivatives on covariant and contravariant vector components are 
defined by 

V i\j = V hj ~ T % v k , (97) 

= w^-H-rj^*, (98) 

where the su mmation on repea t ed in dices is implicit. The symbols 1^ are the Christoffel symbols of the 
second kind 
sections 17.31 and 17.4 



Svnae and SchildL Il978j . Their expressions for the metrics used in this article are given in 



Covariant differentiation of higher order tensors is explained in \Synae and Schildx [19781 ] but we only 
need the rule for a covariant tensor of second order: 

<7ij\k = a i],k - r' fc cry - T jk an . 

If some of the indices of the tensor are contravariant, the rule is changed according to equation (f9"5)) . The 
covariant derivatives of the metric and of the inverse metric are zero: gij\k — and g lJ \ k = 0. 

7.3 Three-dimensional spherical geometry 

The geometry of a thin spherical shell of average radius R can be described with coordinates 9, tp and 
£, respectively representing the colatitude, longitude and radial coordinates. The radial coordinate £ is 
zero on the reference surface (i.e. the sphere of radius R) of the shell. The non-zero components of the 
metric are given by 

g ee = (R + () 2 , 
g vv = (R + C) 2 sin 2 9 , 
■9C C = 1 • 



27 



R + C 



The non-zero Christoffcl symbols are given by 

Ke = -(R + C), 

J< v = -(R + C) sin 2 9, 
7.4 Two-dimensional spherical geometry 

If 9 and <p respectively represent the colatitude and longitude coordinates, the non-zero components of 
the metric on the surface of the sphere are given by 

gee = 1 , 

g w = sin 2 6. (99) 
The non-zero Christoffel symbols are given by 

Y% = -sin 9 cos 9, 

Ko = r e° v = cot0 - 

The double covariant derivatives of a scalar function / are thus given by 

f\e\8 — f,e,e > 

f\e\ip — f\tp\e = f,e,<p — catQf^ , 
f\ip\ip = f,<p,<p + sin 9 cos 9 ■ 
An antisymmetric tensor Sij is defined by 



V det gij e 



n > 



where £y is the antisymmetric symbol invariant und er coordinate transform ations: eg v = —e v g = 1 5 
egg = = (fy is usually called a tensor density; Svnae and Schild 1978| call it a relative tensor of 



weight -1). The non-zero covariant components of are given for the metric of the spherical surface by 

£8ip = -£ v e = sin9 . 

The non-zero contravariant components, = g <7 J are given by 

£ e v = _ £V e = cgc6 ^ 

The covariant derivative of the tensor ey is zero: e^ifc = 0. 
7.5 Gradient, divergence, curl and Laplacian 

Various differential operators on the surface of the sphere can be constructed with covariant derivatives. 
In this section , / and t are scalar functions defined on the sphere and v is a vector tangent to the sphere. 



Backus 



1986] gives more details on surface operators and on Helmholtz's theorem. 
As mentioned in Appendix 17. 2\ the covariant derivative of a scalar function / defined on the sphere 
is a covariant vector tangent to the sphere whose components are f^g and f jlp . The surface gradient of / 
is the same vector with its components expressed in the normalized basis (9, ip): 

V/ = / lfl + csc0/, v 0. (100) 
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The contraction of the covariant derivative with the components of a vector v yields a scalar: 

v\ t = v e e + cot6v 

The surface divergence is the corresponding operation on the vector with its components expressed in the 
normalized basis (0,cp): 

V • v = csc# ((sm6v$) e + v v ,<p^ ■ (101) 

Since the result is a scalar, v % u = V • v. A useful identity is 

V- (/v) = V/-v + /V-v. (102) 

The contraction of the antisymmetric tensor with the covariant derivative of a scalar t yields the 
covariant components of a vector v: 

Vi = 9 jk £ik tj . 

The components are given for the metric (j9"9"|) by vg — csc9t^ and v v = — smdtg. If t is considered 
as the radial component of the radial vector t = t f (the covariant radial component is equal to the 
normalized one), Vi are the non-zero covariant components of the three-dimensional curl of t, which is 
tangent to the sphere. This fact justifies the definition of the surface curl of t, which is equal to the 
vector v but with components given in the normalized basis (0,tp): 

V x t = esc 9 t v - t,g<p. (103) 

The contraction of the double covariant derivative acting on a scalar / defines the surface Laplacian: 

A/ = 9 ij f W 

= f,9M + cot0 l g + esc 2 9f tVtV . (104) 

The surface Laplacian can also be seen as the composition of the surface divergence with the surface 
gradient: A/ = V • V/. 

According to Hclmholtz's theorem, a vector tangent to the sphere can be written as the sum of the 
surface gradient of a scalar / and the surface curl of a radial vector t f : 

v = V/ + Vx(ir). (105) 

While t is alw ays called the toroidal scalar (or potential) for v, there is no standard terminology for / 
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19861 ] calls / the consoidal scalar for v. Some authors [e.g. Banerdt . Il986| ] call / the poloidal 



potential for v. The origin of this use lies in the theory of mantle convection, in which plate tectonics are 
assumed to be driven by mantle flow. Under the assumption of an incompressible mantle fluid, the velocity 
field of the fluid is solenoidal, i.e. its 3-dimensional divergence vanishes. In such a case, the velocity field 
can be decomposed into a po loidal part (V x V x (Ff)) and a toroidal part (V x (Q f )), where differential 



operators are 3-dimensional BackusY Il986l |. If the velocity field i s tangent to the spheric al surface, the 



Forte and Peltier 



19871. However the 



poloidal component at the surface is also the consoidal component 
fields for which we use Helmholtz's theorem, i.e. the tangential surface load and the tangential surface 
displacement, do not belong to 3-dimensional solenoidal vector fields. We thus prefer to use the term 
'consoidal'. 

The surface divergence of v depends only on the consoidal scalar /: 

V-v = A/. (106) 
The two-dimensional version of Gauss theorem is 

dto V • v = . (107) 

s 

where dio = sin# d6 dip and the integral is taken over the whole spherical surface. It can be proven with 
formula (flQljl . 
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7.6 Rigid displacements 



At the surface of a sphere subjected to deformation, the displacement u of a point can be expressed with 
the help of Hclmholtz's theorem (|105|) in terms of three scalar functions (w, S, T) depending on 9 and tp: 



u = wr + VS + V x (Tr) 



(108) 



Strains (and stresses) vanish for rigid displacements. Equations (|69[) show that strains vanish when 
(w, S, T) are of degree one, with S = w (recall that the operators Oi annihilate the degree one). Assuming 
these conditions, we now show that \itransi — w f + V w represents a rigid translation whereas u rot — 
V x (tr) represents a rigid rotation of the sphere. We choose as basis the real spherical harmonics of 
degree one which form the components of the radial unit vector in Cartesian coordinates: 



(109) 



{Y x ,Y y ,Y z ) = (sin 8 cos ip, sin 9 sin tp, cos ( 
= (x,y,z)-?. 



We need the surface gradient of the real spherical harmonics which can be computed with formulas 
and pOO)) : 

(V Y x , V Y y , V Y z ) = (x — sin 9 cos <p f , y — sin 9 sin tp f , z — cos 9 f ) . (HO) 
If the expansion of w in the degree-one basis is w = aY x + bY y + cY z , then 

wf + Vw = flx + !)y + cz, 

so that Utransi is indeed a rigid translation of the sphere. 

If the expansion of T in the degree-one basis is T — a'Y x + b'Y y + c'Y z , then 

V x (T f ) = a' ( — sin tp6 — cos 9 cos <p J + b' I cos tp6 — cos sin <p ) + c' sin <p , 



so that u rot includes a rigid rotation of the sphere, with (a', 6', c') being the angles of rotation around 
the axes (x, y,z), respectively. Though u rot seems to include a uniform radial expansion, one should 
recall that linearized strain-displacement equations are not valid for large displacements. Since strains 
vanish, the radial expan sion is not p hysical and u ro t represents a pure rotation. Finite deformations are 



for example discussed in \LovJi |l944t [pp. 66-73] and I Sokolnikom |l956j [pp. 29-33] 



7.7 Differential identities for the operators Oi 

The differential operators Oi defined by equations (fT5)) satisfy differential identities useful when obtaining 
the flexure equations. They are special cases of differential identities valid in curved spaces. The presence 
of curvature makes the parallel transport of vectors path-dependent; this property quantifies the curvature 
of space and can be expressed as the lack of commutativity of the covariant derivatives of a vector v: 

Vt\j\k - Vi\k\j = Riljk V 1 , (111) 

where Rujk are the covariant components of the Riemann tensor. On the sphere, the Riemann tensor 
has only one independent component that is non-zero, Rg v g v — — sin 2 9. Other components are related 
by the symmetries R a/3l s = -Rp ai S = -Raps-, = R~,8 a p- 

The substitution of f i to Vi in the commutation relation (jllip provides two differential identities 
satisfied by double covariant derivatives acting on scalar functions: 

(cse0f\ v \ v ) e -csc6f\ v \g#-cos0fie\g+8in.6f t g = 0, 
f\e\e,ip — f\<p\o,e — cafc0f\ v \e + f,ip — 0- 
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The replacement in the above equations of the double covariant derivatives by the normalized dif- 
ferential operators (|3"Tj) yields the following identities: 



(sine 2 f) e -(0 3 f) v - cosO o 1 f = (II), (112) 
(sine 3 f) e - (0^)^+ cos e0 3 f = (12). (113) 



These identities can also be directly checked with the definitions (|15p of the operators 0;. 

The identities (I1)-(I2) can be differentiated to generate identities of higher order. A first useful 
identity is obtained from sin 9(11) e — (12) = 0: 

csc 2 9^(sin 2 e(<D 2 f) g ) g + (0 1 f)^-2(sin9 (0 3 /)J ) -cot 9 (Oif) g + 2(D\f = A'/ . (114) 
A second useful identity is obtained from (II) + sin 0(12) g = 0: 
esc 2 9 ((sin 2 e (O a f) t0 ) g - W),^ + (sin 9 ((<% - Oi) /), J J + cot (<%/) _ e - 2 <%/ = . (115) 
7.8 No degree one in operators A and B 



We want to prove that the operators A and B defined by equations (|33|) and (|37|) do not have any 
degree-one term in their spherical harmonic expansion: 



dtu A(a;b)Y* p = (p = -1,0,1), (116) 

J duB{a;b)Y^ = (p = -1,0,1), (117) 

where (a, 6) are arbitrary scalar functions on the sphere, dio = sin 9 d9 dip and the integral is taken over 
the whole spherical surface. 

This property is not a straightforward consequence of constructing A and B with X>,j as a building 
block. Although A and B can be factored into terms without degree one (such as Dy-a or A'a), the 
product of the factors may contain degree-one terms in its spherical harmonic expansion. 

Without loss of generality, we can prove the above identities with the arguments (a, b) being spherical 
harmonics of given degree and order. The general result is then obtained by superposition. Let a and b 
be spherical harmonics of order m and n: a ~ e tm f and b ~ e mv (we will not use their harmonic degree 
in the proof). All derivatives with respect to ip in the operators A and B can then be replaced with the 
rules a j¥ , — > ima and b }V — > in&. The integral over (/? in equations (Ill€i[l - (I117|) gives 

2 \p e i(m+n- P ) V =27r5 m+n _ p , , 

so that the integral is zero unless n = p — m. 

First consider the case p = (n = —m), that is the projection on the zonal spherical harmonic of 
degree one. We thus have to calculate J* d9 Ao and JJ r d9 B$ with 

A = sin 9 cos 9 A(a ; b) , 
£>o = sin9 cos 9 B(a;b) . 

The trick consists in rewriting the integrands as total derivatives: 

A a = (cos 2 ea.fi bfi + cot 6 (sin 2 9 - m 2 ) (ai>),e + (sin 2 9 + m 2 esc 2 6* 00826*) , 

_Bq = — i Tif cos 6* a,e bfi + sin 9 a b t g — esc 9 cos 2 6 1 a_g b + cos 9 a bj. 
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The sought integrals are thus given by 



/' 

Jo 
Jo 



d9 A n = a_gb,g — m cot 6 (ab) t g + m esc 9 ab 
d6 Bq = —im cos 9 a g bg — esc 9 a t g b + cos 8 ab 



where we have dropped the terms containing at least one power of sin# which vanish at the limits; we 
have also replaced cos 2 8 and cos 2$ by their value at the limits. The remaining terms can be evaluated 
by recalling the dependence in sin 8 of the spherical harmonics: a = (sin0)H ao and b = (sin9)' m ' bo, 
where ao and b are polynomials in cos 9. The only non-zero terms at the limits of the integrals are those 
for \m\ = 1, in which case we have at the limits: a,gb g = a bo, cot 8 (ab),g — 2a bo, esc 2 8 ab = ao&o, 
esc 8 a ^gb — cos ao fro, ab : g : g = 0. However these terms cancel in the sums so that the integrals vanish for 
all m. This completes the proof for the case p = 0. 

Now consider the case p = ±1 (n = —m ± 1), that is the projections on the sectoral spherical 
harmonics of degree one. We thus have to calculate f£ dd A±\ and JJ dd B±i with 



B ±1 



skr 8 A(a;b) 



sin 2 8 B {a: 



We again write the integrands as total derivatives: 

A±i — ^ sin cos b g — (sin 6* cos + m 2 ) (ab) g — (cos 2 9 =F 2m) a t g b 
+ sin 2 9ab t g + 2m(m =F 1) cot Sab 



B±i = — i y (m =p 1) sin^a,e b$ + msin# ab t Q t g — (m^ 1) cos# a t $ b 
— m cos 9 ab.g — m(l =p m) esc 9 ab 



The sought integrals are thus given by 



d9A 



±i 



d6B±i = 



m 2 [ab)fi — (1 =F 2m) a^g b + 2m(m =F 1) cot 9 ab 



(m =f 1) cos 9 a y gb + m cos 9 a b t g + m(l ^ m) esc 9 ab 



where we have dropped the terms containing at least one power of sin 9 and replaced cos 2 9 by its value 
at the limits. The remaining terms can be evaluated as in the case p = 0, but with a = (sin6>)' m ' ao and 
b = (sin 0)I" 1= f 1 5 _ All terms give zero at the limits of the integrals for all values of m. This completes 
the proof for the case p = ±1. We have thus proven the identities (|116|) - (|117p . 
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